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Abstract 


This  study  is  aimed  at  development  and  application  of  a  new  3D  wave  propagation  and 
modeling  method  for  regional  waves  in  heterogeneous  crustal  waveguides  using  one-way  wave 
approximation.  The  great  advantages  of  one-way  propagation  methods  are  the  fast  speed  of 
computation,  often  by  several  orders  of  magnitude  faster  than  the  full-wave  finite-difference 
and  finite-element  methods,  and  the  huge  savings  in  internal  memory.  As  the  first  step, 
we  solve  the  2D  SH-wave  problem.  A  half-space  GSP  (Generalized  Screen  Propagator)  is 
formulated  for  the  SH  half  space  problem.  Two  versions  of  the  half-space  GSP  are  derived: 
the  wide-angle  pseudo-screen  and  the  phase-screen.  The  treatment  of  the  Moho  discontinuity 
is  also  discussed.  Comparisons  with  a  wavenumber  integration  method  for  flat  crustal  models 
and  a  finite-difference  algorithm  for  heterogeneous  models  show  excellent  agreement.  For  a 
model  with  propagation  distance  of  250  km,  dominant  frequency  at  0.5  Hz,  the  GSP  method 
is  about  300  times  faster  than  the  finite-difference  method  with  a  similar  accuracy.  These 
comparisons  demonstrate  the  accuracy  and  efficiency  of  the  method.  We  apply  our  method 
to  simulating  regional  wave  propagation  in  different  types  of  complex  crustal  waveguides 
including  those  with  small-scale  random  heterogeneities  and  random  rough  interfaces  of 
sedimentary  layers.  Synthetic  seismograms  and  snapshots  are  shown  to  facilitate  the  study 
of  path  effects  of  Lg  waves.  Slowness  domain  analysis  is  especially  useful  in  investigating 
energy  transfer  in  crustal  waveguides  and  for  determining  which  part  of  the  energy  can  be 
trapped  in  the  waveguide.  The  influence  of  crustal  heterogeneities  and  rough  interfaces  on 
Lg  amplitude  attenuation  and  Lg  coda  formation  are  studied.  It  is  found  that  small  scale 
heterogeneities  and  rough  interfaces  can  scatter  the  Lg  waves  out  of  the  guided  modes  causing 
leakage  into  the  mantle.  The  leakage  loss  due  to  scattering  can  be  an  important  factor  for 
Lg  attenuation  and  blockage. 
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1  Introduction 


The  study  of  path  effects  of  complex  structures  and  heterogeneities  on  the  excitation  and 
propagation  of  regional  phases  in  different  areas  remains  critical  for  many  seismological  prob¬ 
lems,  such  as  the  study  of  crustal  structures  and  the  discrimination  and  yield  estimation  for 
low-yield  nuclear  test  monitoring.  An  ideal  numerical  method  for  investigating  regional  wave 
propagation  should  have  the  capabilities  to  deal  with  crustal  models  that  include  horizontal 
and  vertical  variations  with  scales  from  geological  boundaries  to  small  scatterers,  including 
random  heterogeneities,  as  well  as  handling  3D  Q  structures.  For  the  purpose  of  monitoring 
the  CTBT  (Comprehensive  Test  Ban  Treaty)  at  regional  distances,  simulation  algorithms 
are  desirable  for  generating  synthetic  waveforms  for  high  frequencies  up  to  25  Hz  at  distance 
greater  than  1000  km. 

Substantial  efforts  have  been  made  in  modeling  regional  wave  propagation.  Methods 
based  on  layered  earth  models,  such  as  the  reflectivity  and  mode  summation  methods  (e.g. 
Bouchon  et  al.,  1985,  Kennett,  1989,  1990;  Maupin,  1989;  Baumgardt,  1990;  Campillo, 
1990;  Campillo  and  Paul,  1992;  Campillo  et  al.,  1993;  Gibson  et  al.,  1994)  have  very  high 
efficiency  and  can  be  applied  to  relatively  high  frequencies,  but  they  can  be  used  only  for 
very  simplified  cases  with  layered  or  smoothly  varying  layered  models,  and  lack  the  capability 
to  simulate  more  realistic  crustal  waveguides.  Modeling  techniques  that  can  treat'  realistic 
3D  heterogeneous  media  rather  than  smoothly  varying  layered  media  are  needed  to  test  and 
study  many  observations  and  hypotheses.  Sudden  changes  of  crustal  thickness,  strong  lateral 
variations  and  irregular  3D  heterogeneities  are  among  the  problems  requiring  new  modeling 
methods.  As  pointed  out  by  Campillo  et  al.  (1993),  actual  Lg  amplitudes  were  reduced  more 
than  10  times  for  paths  passing  through  an  anomalous  zone  on  the  east  side  of  the  Alpine 
range,  while  the  modeling  results  using  existing  methods  (including  the  effect  of  known 
large-scale  lateral  structural  variation)  only  account  for  20  -  30  %  of  the  amplitude  reduction. 
Other  mechanisms  such  as  the  scattering  and  attenuation  by  small-scale  heterogeneities  must 
be  taken  into  account. 

Finite-difference  methods  (e.g.,  Xie  and  Lay,  1994;  Goldstein  et  al.,  1996;  Jih  1996)  and 
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pseudo-spectral  methods  (e.g.,  Kosloff  et  al.,  1990;  Orrey  et  al,  1996;  Schatzman,  1996)  are 
very  flexible.  Theoretically  they  can  deal  with  arbitrarily  heterogeneous  media.  However,  the 
capability  of  the  present  day  computers  usually  restricts  them  to  short  propagation  ranges 
and  relatively  low  frequencies,  which  prevents  them  from  being  applied  to  more  realistic 
models. 

Recently,  the  generalized  screen  method  has  been  introduced  into  seismic  wave  simu¬ 
lations.  The  generalized  screen  method  is  based  on  the  one-way  wave  equation  and  the 
one-return  approximation.  The  one-way  wave  propagator  GSP  (Generalized  Screen  Propa¬ 
gator)  neglects  backscattered  waves,  but  correctly  handles  all  the  forward  multiple-scattering 
effects,  e.g.,  focusing/defocusing,  diffraction,  interference,  and  conversion  between  differ¬ 
ent  wave  types.  The  one-return  approximation  is  also  called  the  de  Wolf  approximation 
(De  Wolf,  1971,  1985),  which  neglects  the  reverberation  between  screens  and  can  simulate 
multiple-forescattering-single-backscattering  (MFSB).  Significant  progress  has  been  made  on 
the  development  of  an  Elastic  Complex  Screen  (ECS)  method  for  modeling  elastic  wave  prop¬ 
agation  and  scattering  in  arbitrarily  complicated  structures  (Wu,  1994,  1996;  Xie  and  Wu, 
1995,  1996).  The  method  is  two  to  three  orders  of  magnitude  faster  than  the  elastic  finite- 
difference  method  for  medium  size  3D  problem.  The  screen  method  has  been  successfully 
used  in  forward  modeling  (Wu,  1994;  Xie  and  Wu  1995,  1996;  Wu  and  Huang,  1995)  and 
as  backpropagators  for  seismic  wave  imaging/migration  in  both  acoustic  and  elastic  media 
(e.g.  Wu  and  Xie  1994;  Huang  and  Wu  1996).  In  the  crustal  waveguide  environment,  major 
wave  energy  is  carried  by  forward  propagating  waves,  including  forward  scattered  waves,  and 
therefore  the  neglect  of  backscattered  waves  in  the  propagation  will  not  change  the  main  fea¬ 
tures  of  regional  waves  in  most  cases.  By  neglecting  backscattering  in  the  theory,  the  method 
becomes  a  forward  marching  algorithm  in  which  the  next  step  of  propagation  depends  only 
on  the  present  values  of  the  wavefield  in  a  transverse  cross-section  and  the  heterogeneities 
between  the  two  cross-sections.  The  saving  of  computing  time  and  storage  is  enormous  and 
makes  it  a  viable  method  for  large  3D  elastic  wave  propagation  problems. 

As  the  first  step,  Wu,  Jin  and  Xie  (1997)  introduced  the  GSP  method  to  simulate  and 
investigate  SH  Lg-wave  propagation.  In  this  report,  we  will  discuss  many  applications  of  this 


3 


method  to  complicated  crustal  waveguide  structures.  In  section  2  we  first  give  the  wide-angle 
pseudo-screen  formulas  for  one-way  SH-waves  in  a  half  space.  Then  a  phase-screen  algorithm 
is  obtained  under  the  small- angle  forward  scattering  approximation.  The  treatment  of  the 
Moho  discontinuity  is  also  discussed.  In  section  3,  we  compare  the  simulation  results  of  the 
method  with  those  from  wavenumber  integration  and  finite-difference  methods.  For  a  flat 
crustal  model,  for  which  the  wavenumber  integration  (WI)  method  (reflectivity  method)  is 
considered  as  accurate,  our  result  shows  excellent  agreement  with  the  WI  method.  For  a  het¬ 
erogeneous  crustal  model  with  propagation  distance  of  250  km  and  dominant  frequency  of  0.5 
Hz,  the  GSP  method  is  about  300  times  faster  than  the  finite-difference  method  with  a  sim¬ 
ilar  accuracy.  These  comparisons  demonstrate  the  accuracy  and  efficiency  of  the  half-space 
GSP  method.  In  section  4  we  apply  our  method  to  regional  wave  propagation  in  different 
types  of  complex  crustal  waveguides  including  those  with  small-scale  random  heterogeneities 
and  random  rough  interfaces  of  sedimentary  layers.  The  influence  of  these  heterogeneities 
and  rough  interfaces  on  Lg  amplitude  attenuation  and  Lg  coda  formation  are  shown  to  be 
significant.  Section  5  focuses  on  the  slowness  domain  analysis.  Although  investigating  wave 
propagation  in  space-time  domain  is  more  straightforward,  slowness  domain  analysis  pro¬ 
vides  a  clearer  picture  of  energy  propagating  directions  as  well  as  better  criteria  for  energy 
partitioning.  Conclusion  and  discussions  are  given  in  the  final  section. 

2  Theory 

First  we  derive  the  wide-angle  screen  formula  for  the  2D  SH-wave  one-way  propagation  in 
an  elastic  half  space.  The  equation  of  motion  in  a  linear,  heterogeneous  elastic  medium  can 
be  written  as  (Aki  and  Richards,  1980) 

— u2p(r)u(r)  =  V  •  cr(r)  (1) 

where  u  is  the  frequency,  r  =  (x,  z)  is  a  2D  position  vector,  u  is  the  displacement  vector,  a  is 
the  stress  tensor  (dyadic)  and  p  is  the  density  of  the  medium.  Here  we  assume  no  body  force 
exists  in  the  medium.  For  an  isotropic  2D  elastic  medium,  the  SH  and  the  P-SV  waves  are 
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decoupled.  In  this  article  we  will  treat  only  the  SH  problem  to  demonstrate  the  applicability 
of  screen  propagators  to  crustal  waveguide  problems. 


2.1  General  wide-angle  formulation 

For  a  2D  SH  problem,  only  the  y-component  of  the  displacement  field  exists  and  the  diver¬ 
gence  of  the  stress  tensor  can  be  simplified  into 

V  •  <r(r)  =  V  •  (pVu) 

d  r  ,  .  d  .  d  r  3  ,  ... 

=  &Wr)&ul  +  (2) 

where  u  =  uy,  V  is  the  2D  gradient  operator,  and  p  is  the  shear  rigidity.  We  decompose  the 
parameters  of  the  elastic  medium  and  the  total  wave  field  into 

p  =  po  +  Sp 
p  =  p0  +  Sp 

u  =  u°  +  U  (3) 

where  p0  and  po  are  the  parameters  of  the  background  medium,  Sp  and  Sp  are  the  corre¬ 
sponding  perturbations,  u°  is  the  primary  field  and  U  is  the  scattered  field.  Then  the  SH 
wave  equation  can  be  rewritten  as 

PqS72u  +  u>2p0u  =  —  [u25pu  +  V  •  Sp'Vu]  ,  (4) 


or 

(V2  +  k2)u{ r)  =  —  k2F(r)u(r)  ,  (5) 

where  k  =  uj/v  is  the  wavenumber  in  the  background  medium  and  v  is  the  background  S 
wave  velocity  defined  by 

v  =  yj  Po/ Po  (6) 

In  the  right  hand  side  of  (5),  F(r)  is  a  perturbation  operator 


with 


£p(r)  = 

£p{r)  = 


6p{r) 

Po 

Mr) 

Po 


(8) 

(9) 


Equation  (5)  is  a  scalar  Helmholtz  equation.  With  a  half-space  scalar  Green’s  function 
gh,  the  scattered  field  U  can  be  written  as 


U(ri)  =  k2  f  d2Tgh(r1-,r)F(r)u(T)  ,  (10) 

J  V 

where  the  2D  volume  integration  is  over  the  volume  V  including  all  the  heterogeneities  in 
the  modeling  space.  Using  the  Gauss  divergence  theorem  the  above  equation  reduces  to  (cf. 
Wu,  1994) 

Ufa)  =  k2  f  d2r  {/(ri;r)ep(r)n(r) 

J  V  V 

-  ^V/(rx;r) -eM(r)Vw(r)J  (11) 

Under  the  forward-scattering  approximation,  or  more  generally  the  multiple-forescattering- 
single-backscattering  (MFSB)  approximation  (De  Wolf,  1971,  1985;  Wu  and  Huang,  1995; 
Wu,  1996),  the  total  field  and  Green’s  function  under  the  integration  in  the  above  equation 
can  be  replaced  by  their  forward-scattering  approximated  counterparts,  and  the  field  can 
be  calculated  by  a  one-way  marching  algorithm  along  the  x-direction  using  a  dual  domain 
technique.  We  will  derive  in  the  following  the  dual  domain  expressions  of  the  marching 
algorithm  for  the  heterogeneous  half-space  case. 

For  each  step  of  the  marching  algorithm  under  the  forward-scattering  approximation,  the 
total  field  at  oq  is  calculated  as  the  sum  of  the  primary  field  which  is  the  field  free-propagated 
in  the  half-space  from  x'  to  aq,  and  the  scattered  field  caused  by  the  heterogeneities  in  the 
thin-slab  between  x'  and  aq.  The  thickness  of  the  slab  should  be  made  thin  enough  to  ensure 
the  validity  of  the  local  Born  approximation.  Under  this  condition,  the  Green’s  function 
can  be  approximated  by  the  homogeneous  half-space  Green’s  function.  The  latter  can  be 
obtained  by  the  image  method  (Morse  and  Feshbach,  1953).  The  stress  should  vanish  at  the 


6 


free  surface  z  =  0,  leading  to  jxdu/dz  =  0,  a  Neumann  boundary  condition.  Therefore,  we 
have 

Po(ri;r)  =  5o(ri;r)  +p0(ri;r*)  (12) 


where  g0  is  the  infinite  homogeneous  Green’s  function  and  r*  is  the  mirror  point  of  r  with 
respect  to  the  free  surface. 

Taking  the  Fourier  transform  of  Eq.  (11)  along  Zi,  for  the  case  of  a  thin-slab  we  have 


U(xuKz) 


rx i  roo 

=  k2  dx  dz{g^ (aq,  Kz\  x,  z)ep(x,  z)u0(x,  z) 
J  xf  J  0 

-  z)  ■  £„(*,  z)Vu(x,  z)} 


(13) 


We  know  the  free  space  Green’s  function  in  wavenumber  domain  (Wu  and  Huang,  1995;  Wu, 
1996) 

go(x1,Kz;x,z)  =  i-x)e-iK,z  (14) 

27 

with 

7  =  -  Kl  ,  (15) 

Therefore, 


go(xi,Kz-,x,z)  =  go(xi,Kz;x,z)  +  g0{x1,Kz-,x,  -z) 

—  _Le*7(zi-z)je ~iKzz  e+iKzzj 
27 


=  JLei7(*i-x)2c0S  (KZZ) 

27 


In  a  similar  way  we  can  obtain 


and 


Vg0{xi,Kz-,x,z)  =  ^-(jex  +  Kzez)elj(xi  x)e  tKzZ 
27 


+ 

=  iet7(zi_x){e22cos (Kzz)  -  ez{Kz/j)2ism(Kzz)} 


where  ex  and  ez  are  unit  vectors  in  the  x  and  2  directions,  respectively. 


(16) 


(17) 


(18) 


7 


Substitute  g[}  and  V#q  into  (13),  leading  to 


U{xuKs) 
U p(x  i,  Kz) 

Up(x  i,  Kz) 


Up{xuKz)  +  Up{xi,Kz) 

jU  rx  i  r  oo  Zr* 

—  /  dxell{-Xl~x^  /  dz2  cos^z) -ep(x,  z)u0(x,  z) 
2  Ji'  Jo  7 


-  —  P  dxe^Xl~x) 

2  Jx' 

r°°  -  K  - 

/  d;z[2cos(/if*z)awo  —  2isin(/srz2:) — ^9zuo]e„(rr,  z) 
Jo  7 


where 


«  _  1  a_ 

1  <9rr 

J__d 

ik  dz 


a  = 


are  dimensionless  partial  derivatives.  Or  we  can  write  these  equations  as 

E/,(*i,K,)  =  £  [Z'  dxe^“-‘>C{-ef(z)Mz)} 
z  Jx'  7 

ik  rx  1  f  _  _ 

Kz)  =  ~—  /  <£rei7(xi_l)  Cfe^auofc)]  -  *5[— e^(«)a«o(^)] 
z  Jx'  f  7  . 

where  C[f(z)]  and  S[f(z)]  are  cosine  and  sine  transforms,  defined  by 


r  00 

C[f(z)}  =  /  dz2cos(Kzz)f(z) 

Jo 

f  00 

S[f(z )]  =  /  dz2sm(Kzz)f(z) 

Jo 

for  forward  transforms,  and 

C~l[f(Kz)}  =  —  /  dKz2cos(Kzz)f(Kz) 

Z7T  «/  0 
1  f00 

■S_1[/(^)l  =  IT  dK,2sm(Kzz)f(Kz) 

I/K  J  0 

for  inverse  transforms. 

In  Eqs.  (21)  and  (22),  w0,  a^o  and  a«0  can  be  calculated  by 

u0(x,z)  =  J- [°°  a(x',K'z) 

Z7T  J —00 


(19) 


(20) 


(21) 

(22) 


(23) 


(24) 


(25) 
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and 


4«o  (x,z)  = 

8,«o(*,z)  =  iS-'le^-'&u o(x',K',)\  (26) 

Eqs.  (21),  (22),  (25)  and  (26)  are  the  dual-domain  expressions  of  the  wide-angle  screen 
propagator  for  half-space  SH  problems. 

2.2  Implementation  procedure  for  the  half-space  wide-angle  screen 
propagator 

Under  the  forward-scattering  approximation  we  can  update  the  total  field  with  a  marching 
algorithm  in  the  forward  direction.  The  half-space  model  can  be  sliced  into  thin-slabs  per¬ 
pendicular  to  the  propagation  direction.  Weak  scattering  condition  holds  for  each  thin-slab. 
For  each  slab-step  forward,  the  forward-scattered  field  by  the  thin-slab  is  calculated  and 
added  to  the  primary  field  so  that  the  updated  field  becomes  the  incident  field  for  the  next 
thin-slab.  The  procedure  can  be  summarized  as  follows. 

1.  Cosine  transform  the  incident  fields  at  the  entrance  of  each  thin-slab  into  wavenumber 
domain. 

2.  Free  propagate  in  wavenumber  domain  and  calculate  the  primary  field  and  its  gradient 
within  the  slab. 

3.  At  each  horizontal  position  within  the  slab,  inverse  cosine/sine  transform  the  primary 
field  and  its  gradients  into  space  domain  (25  and  26),  and  then  interact  with  the  medium 
perturbations  ep  and  ep. 

4.  Cosine/sine  transform  the  distorted  fields  into  wavenumber  domain  and  perform  the 
divergence  operations  to  get  the  scattered  fields  (21  and  22). 

5.  Calculate  the  primary  field  at  the  slab  exit  and  add  to  the  scattered  field  to  form  the 
total  field  as  the  incident  field  at  the  entrance  of  the  next  thin-slab. 

6.  Continue  the  procedure  iteratively. 
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2.3  Small  angle  approximation  and  the  phase-screen  propagator 

When  the  energy  of  crustal  guided  waves  is  carried  mainly  by  small-angle  waves  (with  respect 
to  the  horizontal  direction),  the  small  angle  approximations  can  be  invoked  to  simplify  the 
theory  and  calculations.  Let  us  first  consider  the  calculations  for  U 'p.  Substitute  (26)  into 
(22),  resulting  in 

l MxuKz)  =  ~  r  f°°  dK'^'^r  dz 

2  J xr  Ztv  Jo  Jo 

^2cos(Kzz)ep{z)2cos(K'zz)(j^) 

+  2sm(Kzz)£tl(z)2sin(K'zz)(^-)(~)^uo(x',K'z )  (27) 

This  is  the  local  Born  scattering  in  wavenumber  domain.  The  above  equation  can  be  written 
as 

U^x i,  Kz)  =  A^uoix',  K'z)  (28) 

where  is  an  operator.  In  the  discrete  form,  it  is  a  matrix,  whose  elements  are  incident 
and  outgoing  (scattered)  wavenumbers  dependent.  For  small-angle  waves,  Kz  <C  7  «  7'  «  k. 
The  second  term  of  Eq.  (27)  is  a  second  order  small  quantity  compared  with  the  first  term, 
and  therefore  can  be  neglected  in  this  case.  Then  we  have 

U(xi,Kz)  =  Up(xi,  Kz)  +  Ull(x1,Kz) 

jlc  .  ,  1  roo  r 00 

=  —  e’7(Xl-x /  dK'z 2  cos(K'zz)  /  dz2cos(Kzz) 

2  2tt  J  0  J 0 

Jx,  d®ei(7"7,)(*-x')  |(*)ep(s)  -  (jM*)  u0(x',K'z) 

«  ikJ^Xl-x,)C[Ss(z)u0(x',z)\  (29) 


where 


1  rXi 

Ss(z )  =  -/  dx[ep(x,z) -e^faz)] 

Z  Jx’ 

~  Arre^) 


(30) 


where  es{z)  is  the  average  S-wave  slowness  perturbation  over  the  thin-slab  at  depth  z, 
£s{z)  =  ^7  lx'1  dx^J--30-  with  s(x,z)  =  5^,  Ax  =  (xi  -  x')  is  the  thin-slab  thick¬ 
ness.  Here  Ss  is  the  corresponding  slowness  screen  for  the  half-space  thin-slab.  Eq.  (29) 
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is  the  screen  approximation  of  the  half-space  SH  problem.  Summing  up  the  primary  and 
scattered  fields  and  invoking  the  Rytov  transform  result  in  the  dual-domain  expression  of 
phase-screen  propagator  for  this  case 

u(xi,Kz )  =  u0(xuKz)  +  U(xl,Kz) 

_  et7(zi— z ')  f  dz2cos(Kzz)[l  +  ikSs(z)]uo(x',z) 

Jo 

«  e**x'~^C  [eikS^z)u0(x' ,z)]  (31) 

where  elkSa ^  is  the  phase  delay  operator. 

2.4  Implementation  procedure  for  the  half  space  phase-screen 
propagator 

Under  the  phase-screen  approximation,  the  heterogeneous  half-space  is  represented  by  a 
series  of  half-screens  embedded  in  the  homogeneous  background  half-space.  The  wave  prop¬ 
agates  between  screens  in  the  wavenumber  domain  and  interacts  with  the  phase-screens  in 
the  space  domain.  The  interaction  is  only  a  phase-delay  operator  (multiplication  in  space 
domain).  The  procedure  can  be  summarized  as  follows. 

1.  Cosine  transform  the  incident  field  at  the  starting  plane  into  wavenumber  domain  and 
free  propagate  to  the  screen. 

2.  Inverse  cosine  transform  the  incident  field  into  space  domain  and  interact  with  the 
shear  slowness  screen  (phase-screen)  to  get  the  transmitted  field. 

3.  Cosine  transform  the  transmitted  field  into  wavenumber  domain  and  free-propagate 
to  the  next  screen. 

4.  Repeat  the  propagation  and  interaction  screen-by-screen  to  the  boundary  of  the  model 
space. 

2.5  Treatment  of  the  Moho  discontinuity 

The  Moho  discontinuity  can  be  treated  in  two  ways.  One  is  to  put  the  impedance  boundary 
conditions  in  the  formulation,  the  other  is  to  treat  the  parameter  changes  as  perturbations 
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Figure  1:  Comparison  of  synthetic  seismograms  along  the  surface  calculated  by  the  screen 
method  (thick  lines)  and  reflectivity  method  (thin  lines)  for  a  flat  crustal  model  (32  km 
thick).  The  source  function  is  a  Ricker  wavelet  with  dominant  frequency  of  1.0  Hz. 

and  therefore  be  incorporated  into  the  screen  interaction.  The  former  has  the  advantage  of 
computational  efficiency.  The  latter  has  the  flexibility  of  handling  irregular  interfaces.  In 
this  paper,  we  adopt  the  latter  approach  and  check  the  validity  of  perturbation  approach  for 
the  Moho  discontinuity  by  a  finite-difference  algorithm. 
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3  Numerical  verification  and  tests 


To  test  the  validity  of  the  half-space  GSP  method  (or  simply  call  it  the  screen  method),  we 
have  conducted  extensive  numerical  tests  and  checked  the  results  with  various  well  known 
numerical  methods,  such  as  the  wavenumber  integration  and  finite-difference  methods. 

First,  in  Figure  1  we  show  the  accuracy  of  the  method  by  comparing  the  synthetic  seis¬ 
mograms  generated  by  the  screen  method  (thick  lines)  with  those  calculated  by  a  reflectivity 
method  (thin  lines)  for  a  flat  crustal  model.  The  crust  has  a  thickness  of  32  km  and  a 
shear  wave  velocity  of  3.5  km/s.  The  mantle  beneath  the  crust  has  a  shear  velocity  of  4.5 
km/s.  The  source  function  is  a  Ricker  wavelet  with  a  dominant  frequency  of  1.0  Hz.  Our 
results  agree  very  well  with  the  reflectivity  method  which  is  considered  very  accurate  for  flat 
layered  media.  The  only  exceptions  are  for  near  vertical  reflections  at  very  small  epicentral 
distances  where  the  screen  method  has  low  accuracy  for  extremely  large  scattering  angles 
with  respect  to  the  propagation  direction.  However,  since  regional  seismograms  are  usually 
recorded  at  large  distances,  this  limitation  does  not  pose  any  real  problem  for  its  application. 
Next,  we  show  the  accuracy  of  the  method  by  comparing  synthetic  seismograms  generated 
by  this  method  with  those  generated  by  a  finite-difference  algorithm  (Xie  and  Lay,  1994). 
For  the  finite-difference  method,  a  fourth-order  elastic  SH-wave  code  is  used  to  calculate 
the  synthetic  seismograms.  The  spatial  sampling  interval  is  0.125  km  in  both  vertical  and 
horizontal  directions,  and  the  time  interval  is  0.015  second.  For  the  screen  method,  the 
spatial  sampling  interval  is  0.25  km  in  the  vertical  direction  and  the  screen  interval  is  1.0 
km.  A  Gaussian  derivative  is  used  as  the  source  time  function  for  both  methods.  Because  of 
the  computational  intensity  of  the  finite-difference  method,  we  did  the  comparison  at  short 
propagation  distances  (250  km)  and  a  relatively  low  frequency  (/0  ~  0.5  Hz).  Figure  2 
gives  the  comparison  between  synthetic  seismograms  from  the  screen  method  and  from  the 
finite-difference  method  for  the  same  flat  layered  model.  Their  agreement  is  excellent. 

Next  we  show  the  comparison  for  a  heterogeneous  crustal  model.  On  the  top  of  Figure  3  is 
the  crust  model  with  a  narrow  passage  (‘neck’  type)  used  to  calculate  synthetic  seismograms. 
On  the  bottom  are  the  synthetic  seismograms  along  a  vertical  profile  at  an  epicentral  distance 
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Figure  2:  Comparison  of  synthetic  seismograms  along  the  surface  calculated  by  the  screen 
method  (thick  lines)  and  the  fourth-order  finite-difference  method  (thin  lines)  for  a  flat 
crustal  model  (32  km  thick).  The  source  function  is  a  Gaussian  derivative  (/0  ~  0.5  Hz). 
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THINNING  CRUSTAL  MODEL 


Vertical  Distance  (km) 


SP(thick  solid  line)  &  FD(thin  line)  for  inhomogeneous  model 


Figure  3:  Comparison  of  synthetic  seismograms  along  a  vertical  profile  at  the  distance  of 
250  km  calculated  by  the  screen  method  (thick  lines)  and  a  finite-difference  method  (thin 
lines)  for  a  laterally  varying  crustal  model  shown  on  the  top  panel. 
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of  250  km.  The  thin  lines  are  from  the  finite-difference  method  and  the  thick  lines  are  from 
the  screen  method.  The  source  is  located  at  a  depth  of  2  km.  Excellent  agreement  can  be 
seen.  For  this  example,  the  GSP  method  took  20  CPU  minutes  on  a  SPARC-4.  To  obtain 
a  similar  accuracy,  the  finite-difference  calculation  took  3000  CPU  minutes  on  a  SPARC-20. 
The  speed  factor  is  about  300. 

We  have  performed  an  accuracy  comparison  between  the  screen  method  and  finite- 
difference  method  for  different  grid  spacings.  It  was  found  that  when  the  grid  spacing 
for  the  finite-difference  method  is  0.25  km,  although  it  still  satisfies  the  stability  criterion 
for  the  finite-difference  method,  there  are  considerable  discrepancies  between  the  results  of 
the  two  methods  (see  the  figure  in  Wu  et  al.,  1997).  This  is  due  to  the  numerical  disper¬ 
sion  of  the  finite-difference  method.  With  a  finer  grid  spacing  of  0.125  km,  the  results  from 
finite-difference  converged  to  those  of  the  GSP  method.  This  shows  that  for  wave  propaga¬ 
tion  through  a  long  crustal  waveguide,  the  finite-difference  method  requires  finer  grid  than 
the  conventional  stability  criterion.  The  grid  spacing  (Az)  of  the  screen  method  for  the 
comparison  is  0.5  km,  and  the  screen  interval  (Ax)  is  1  km. 

4  Investigating  Lg-wave  propagation  in  complex  crustal 
structures 

In  this  section  we  show  some  examples  demonstrating  the  potential  of  the  method  applied 
to  various  problems  of  regional  wave  propagation. 

4.1  Long  range  Lg  propagation  simulation  in  complex  crustal  waveg¬ 
uides 

Visualize  the  path  effects  by  snap  shots.  The  high  efficiency  of  the  method  permits 
the  simulation  of  Lg  wave  propagation  through  long  crustal  wave  guides.  The  snap-shots 
can  be  easily  generated  to  visually  analyze  the  path  effects  of  Lg  wave  propagation.  Figures 
4,  5  and  6  show  the  snap  shots  calculated  by  the  screen  method  for  different  waveguides 
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Comparison  of  Wave  Propagation 
in  Various  Crustai  Waveguides 

( t  =  30  sec) 


Figure  4:  Snap  shots  at  30  sec.  The  development  of  mantle  waves  and  head  waves  can  be 
seen  clearly. 
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Con pari son  of  Uaue  Propagation 
in  Uarious  Crustal  Uaue  Guides 


(  t  =  50  sec) 


Figure  5:  Snap  shots  at  50  sec.  The  development  of  mantle  waves  and  head  waves,  and  the 
influence  of  crust  thining  and  thickening  can  be  seen  clearly. 
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Comparison  of  Wave  Propagation 
in  Various  Crustal  Waveguides 

( t  =  70  sec ) 


Figure  6:  Snap  shots  at  70  sec.  The  development  of  mantle  waves  and  head  waves,  and  the 
influence  of  crust  thining  and  thickening  can  be  seen  clearly. 
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at  30,  50  and  70  sec.,  respectively.  In  each  figure,  from  top  to  bottom  are  results  for  flat, 
narrowing  and  broadening  crustal  waveguides  respectively.  The  source  is  located  on  the  left 
boundary  at  a  depth  of  2  km.  The  development  of  mantle  wave  and  head  wave,  as  well  as 
the  formation  of  crustal  guided  waves  as  multiple  reflections  between  free  surface  and  Moho 
discontinuity  can  be  clearly  seen.  For  the  inhomogeneous  models,  wave  diffraction,  leakage 
to  the  mantle,  wavefront  distortion  and  increase  of  wavefield  complexity  can  also  be  seen 
clearly.  From  the  comparison  it  is  seen  that  the  passage  of  a  narrow  crustal  waveguide  (‘neck’ 
type)  has  greater  effect  on  Lg  leakage  than  the  broad  passage  (‘belly’  type).  In  the  latter 
case  although  the  wavefronts  are  complicated  due  to  scattering  at  the  edges,  there  is  more 
energy  trapped  in  the  crust  than  the  case  of  narrow  passage  in  which  a  large  percentage  of 
energy  leaks  into  the  mantle.  This  example  demonstrates  the  potential  of  the  method  as  a 
tool  for  investigating  the  path  effects  of  different  crustal  structures. 

Long-range  high-frequency  synthetic  seismograms.  The  following  example  shows 
the  capability  of  this  method  for  long  range  high-frequency  synthetic  seismograms  in  a  lat¬ 
erally  varying  structure.  Figure  7  shows  the  laterally  varying  crust  model  used  in  the  calcu¬ 
lation.  Figure  8  shows  the  high-frequency  synthetic  seismograms  on  the  surface  at  distances 
up  to  1000  km.  The  center  frequency  fc  =  5  Hz  with  the  maximum  frequency  fmax  =  10 
Hz.  In  comparison,  the  low-frequency  (/c  =  1  Hz,  fmax  =  2  Hz)  synthetic  seismograms  are 
shown  in  Figure  9.  It  is  clear  that  without  high-frequency  content,  many  of  the  distinctive 
features  associated  with  Lg  measurements  can  not  be  adequately  modeled.  In  other  words, 
a  proper  simulation  method  with  the  capability  to  generate  accurate  high-frequency  signals 
is  a  necessity  for  the  purpose  of  investigating  regional  phases. 

4.2  The  influences  of  random  heterogeneities  and  rough  interfaces 

The  importance  of  small-scale  random  heterogeneities  to  seismic  wave  propagation  is  well 
known.  There  are  extensive  publications  on  this  subject  in  seismology.  However,  the  role 
of  random  heterogeneities  in  Lg  excitation,  propagation,  attenuation  and  blockage,  is  still 
unclear  due  to  the  complexity  of  the  problem.  The  theory  of  wave  propagation  in  unbounded 
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Parameters  of  Crustal  Model 


Layer 

Vs(km/sec) 

Density(g/cm ) 

Thickness(km) 

1 

3.00 

2.60 

1.00 

2 

3.46 

2.80 

14.00 

3 

3.76 

3.00 

22.00 

4 

4.65 

3.30 

Half-Space 

Crustal  Model 


Figure  7:  An  inhomogeneous  crustal  model  used  in  the  calculation  of  high-frequency  syn¬ 
thetic  seismograms.  Shown  in  the  upper  panel  are  model  parameters  and  the  lower  panel 
gives  the  geometry  of  the  model.  The  receivers  are  on  the  surface  and  shown  by  triangles. 
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Synthetic  seismograms  (fc«5,OHz.fmax=  1 0.OHz),  Time  (Sec) 
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Figure  8:  High-frequency  (/c  =  5  Hz,  fmax  =  10  Hz)  synthetic  seismograms  on  the  surface 
at  distances  up  to  1000  km  for  an  inhomogeneous  crustal  waveguide 


Synthetic  seismograms  (fc— 1  .OHz.lmax~2.OHz),  Tlme(Sec) 


Figure  9:  Low-frequency  (/c  =  1  Hz,  fmax  =  2  Hz)  synthetic  seismograms  on  the  surface  at 
distances  up  to  1000  km  for  an  inhomogeneous  crustal  waveguide 
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random  media  has  been  well  developed.  However,  for  waves  in  complex  crustal  waveguides 
with  random  heterogeneities,  the  theoretical  difficulties  are  overwhelming,  and  no  analytical 
tools  are  available  for  performing  realistic  calculations.  Numerical  simulation  is  an  attrac¬ 
tive  alternative  to  the  theory.  Some  finite-difference  simulations  have  been  conducted  (e.g. 
Frankel  and  Clayton,  1986;  Frankel,  1989;  Xie  and  Lay,  1994;  Jih,  1996).  Limited  by  the 
computation  power,  however,  the  finite-difference  results  are  often  for  short  distances  or 
low  frequencies.  Liu  and  Wu  (1994)  have  done  some  numerical  simulation  using  the  phase- 
screen  method,  but  the  media  simulated  are  limited  to  unbounded  media.  The  development 
of  the  half-space  GSP  method,  enables  us  to  simulate  high-frequency  waves  propagating  in 
complex  crustal  waveguides  to  long  distances.  In  the  following,  we  will  show  two  examples 
demonstrating  the  capability  of  the  method. 

Figure  10  shows  a  heterogeneous  crustal  model  representing  a  ‘mountain  root’  with  small- 
scale  random  heterogeneities.  On  the  top  panel  is  the  velocity  model,  and  the  comparisons 
between  synthetic  seismograms  with  and  without  random  heterogeneities  are  shown  on  the 
middle  and  bottom  panels,  respectively.  The  heterogeneities  have  an  exponential  correlation 
function,  with  the  scale  length  ax  =  az  =  1.6  km  (in  horizontal  and  vertical  directions, 
respectively).  The  RMS  velocity  perturbation  is  5%.  The  dominant  frequency  of  the  source 
function  is  2  Hz.  Figures  11  A  and  B  show  the  comparison  between  snapshots  for  waves 
passing  through  the  ‘mountain  root’  with  and  without  random  heterogeneities,  respectively. 
We  see  that  random  heterogeneities  drastically  increase  the  leakage  of  waves  to  the  mantle 
and  the  complexity  of  the  waveforms.  Extensive  numerical  experiments  will  be  conducted  to 
study  the  different  influences  of  various  kinds  of  random  heterogeneities.  It  has  been  shown 
that  the  crustal  random  heterogeneities  are  highly  anisotropic  in  scale  length  (Levander  and 
Holliger,  1992;  Holliger  and  Levander,  1992;  Wu  et  ah,  1994).  The  influences  of  the  random 
heterogeneities  with  different  stochastic  characteristics  will  be  explored  systematically. 

Figure  12  shows  the  comparison  of  synthetic  seismograms  and  snapshots  for  models  with 
and  without  a  rough  interface.  Shown  on  the  first  panel  is  a  crustal  model  with  a  1  km  thick 
low-velocity  top  layer.  The  bottom  of  the  low-velocity  layer  is  a  rough  interface  with  0.2  km 
RMS  random  depth  fluctuations.  The  randomness  has  an  exponential  correlation  function 
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Time(s) 

Seismogram  without  random  media 


Figure  10:  A  heterogeneous  crustal  model  representing  a  mountain  root  with  small-scale  ran¬ 
dom  heterogeneities  (top  panel).  The  comparison  between  synthetic  seismograms  with  and 
without  random  heterogeneities  are  shown  on  the  middle  and  bottom  panels,  respectively. 
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A.  Snapshot  for  waves  passing  a  mountain  root  without  random  heterogeneities 
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B.  Snapshot  for  waves  passing  a  mountain  root  with  random  heterogeneities 

Figure  11:  Comparison  between  snapshots  for  waves  passing  through  a  ‘mountain  root’  with 
or  without  random  heterogeneities,  shown  on  A  and  B,  respectively. 
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and  a  horizontal  scale  length  of  0.5km.  The  second  and  third  panels  show  the  comparison 
of  synthetic  seismograms;  the  bottom  two  panels  show  the  comparison  of  snapshots.  We 
see  that  rough  interfaces  of  sedimentary  layers  can  also  increase  the  mantle  leakage  and 
waveform  complexity  of  regional  waves. 

5  Energy  transfer  and  partitioning 

While  the  wave-field  snapshots  discussed  previously  are  very  intuitive  for  assessing  general 
feature  of  Lg-wave  propagation,  it  is  rather  difficult,  in  a  complex  model,  to  investigate  the 
energy  transport  by  directly  analyzing  waveforms  and  snapshots.  In  this  section,  we  will 
focus  our  investigation  on  wave  energy  transfer  and  partitioning  in  heterogeneous  crustal 
waveguides. 

The  upper  boundary  of  the  crustal  wave  guide  is  the  free  surface,  which  is  a  perfect 
reflector.  The  lower  boundary  of  the  wave  guide  is  the  Moho  discontinuity.  For  waves 
incident  on  the  Moho  discontinuity,  part  of  the  energy  will  leak  into  the  upper  mantle. 
However,  for  waves  incident  on  the  Moho  with  post-critical  angles,  total  reflections  occur 
and  all  the  energy  is  reflected  and  trapped  in  the  waveguide.  Generally  speaking,  the  guided 
energy 

Eg=[  \u(Kz)\2dKz  (32) 

where  Kc  is  the  critical  wavenumber.  Scattering  processes  can  redistribute  the  energy  in 
wavenumber  domain,  causing  the  trapped  energy  leak  to  the  upper  mantle.  Therefore  in¬ 
vestigating  energy  distribution  versus  the  propagation  angle,  or  equivalently  their  vertical 
slowness,  can  give  a  clear  picture  of  which  part  of  the  energy  can  be  trapped  in  the  wave 
guide  and  which  part  of  the  energy  will  leak  into  the  upper  mantle.  We  will  use  a  set  of 
numerical  examples  and  slowness  analysis  to  investigate  the  interactions  between  Lg-waves 
and  different  types  of  wave  guides. 
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Snapshot  at  T=80s  with  rough  interface 


Figure  12:  Influence  of  a  rough  interface  on  Lg  propagation.  The  top  panel  is  a  crustal 


model  with  a  sedimentary  layer  for  which  the  bottom  is  a  rough  interface.  The  second  and 


third  panels  show  the  comparison  of  synthetic  seismograms;  the  bottom  two  panels  show  the 


comparison  of  snapshots. 
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Angular  Spectrum  for  crustal  waveguides  at  f=[0.6,1.9] 


Distance  (km) 

100  200  300  400  500  600 


Figure  13:  Energy  distribution  for  a  uniform  crust  model.  From  top  to  bottom,  crustal 
structure;  energy  distribution  versus  vertical  slowness  and  distance;  and  energy  attenuation 
versus  distance. 
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Angular  Spectrum  for  crustal  waveguides  at  f=2.0Hz 
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Angular  Spectrum  for  crustal  waveguides  at  f=1 ,5Hz 
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Angular  Spectrum  for  crustal  waveguides  at  f=1  .OHz 
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Angular  Spectrum  for  crustal  waveguides  at  f=0.5Hz 


Figure  14:  Energy  distribution  of  a  uniform  crust  model  versus  vertical  slowness  and  distance 
for  different  frequencies.  From  top  to  bottom,  the  frequencies  are  2.0,  1.5,  1.0  and  0.5  Hz 
respectively. 


5.1  Deterministic  waveguide  models 

Figure  13  shows  the  result  for  a  simple  waveguide.  The  model  is  composed  of  a  30  km 
thick,  horizontal,  constant  velocity  crust  over  a  constant  velocity  upper  mantle.  The  S-wave 
velocities  and  densities  are  V\=  3.5  km/s,  pi  =  2.7  g/crr?  and  V2=  4.5  km/s,  p2—  3.3  g/cm3 
for  the  crust  and  upper  mantle  respectively.  Hereafter,  we  will  refer  to  this  model  as  the 
reference  model.  Shown  in  the  upper  panel  of  Figure  13  is  the  waveguide  structure.  In 
the  middle  is  the  energy  distribution  versus  the  vertical  slowness  and  distance  up  to  600 
km.  The  vertical  coordinate  is  the  normalized  vertical  slowness  Kz/k,  corresponding  to  the 
cosine  of  incident  angles.  Note  that  zero  vertical  slowness  means  horizontal  propagation. 
The  frequency  range  is  from  0.6  to  1.9  Hz.  At  the  initial  stage,  there  is  considerable  energy 
with  large  vertical  slowness,  i.e.,  with  steep  angles.  After  multiple  reflections,  energy  with 
larger  vertical  slowness  is  depleted  due  to  the  leakage  to  the  mantle,  leaving  the  energy  with 
small  vertical  slowness,  i.e.,  the  guided  waves,  propagating  in  the  waveguide.  Shown  in  the 
lower  panel  is  the  wave  energy  versus  the  distance.  The  energy  is  calculated  from  synthetic 
seismograms  on  the  free  surface.  It  can  be  compared  with  the  observed  Lg  waves  from 
seismic  networks.  After  the  first  100  km,  the  energy  is  basically  constant.  Note  that  this 
result  is  for  a  two-dimensional  range  independent  waveguide  without  intrinsic  or  scattering 
attenuation.  For  a  real  3D  model,  the  geometrical  spreading  factor  should  be  different  from 
here.  Figure  14  shows  the  energy  distribution  versus  the  vertical  slowness  and  distance  for 
the  reference  model  for  individual  frequency  components.  From  top  to  bottom  are  plots 
for  frequencies  2.0,  1.5,  1.0  and  0.5  Hz,  respectively.  The  energy  distribution  shows  some 
patterns  in  the  slowness-distance  domain.  These  result  from  the  interference  of  waves  from 
the  source  as  well  as  from  the  free  surface  and  bottom  reflections.  After  a  certain  distance, 
the  patterns  are  relatively  stable,  that  is,  certain  modes  have  formed  in  the  waveguide  and 
they  carry  the  energy  to  a  long  distance.  Later  we  will  see  the  waveguide  structures,  in 
both  vertical  and  horizontal  directions,  can  affect  the  propagating  mode,  i.e.,  the  patterns 
of  energy  distributions. 

The  second  example,  shown  in  Figure  15,  is  similar  to  the  first  one  except  there  is  a  crust 
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Figure  15:  Energy  distribution  for  model  with  a  crustal  necking.  From  top  to  bottom,  crustal 
structure;  energy  distribution  versus  vertical  slowness  and  distance;  energy  attenuation  ver¬ 
sus  distance. 
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Angular  Spectrum  for  bellying  crustal  model 
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Energy  attenuation  versus  distance  (solid  line-bellying) 


Figure  16:  Energy  distribution  for  model  with  a  crustal  thickening.  From  top  to  bottom, 
crustal  structure;  energy  distribution  versus  vertical  slowness  and  distance;  and  energy  at¬ 
tenuation  versus  distance. 
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necking  in  the  waveguide  where  the  crust  thickness  decreased  from  30  km  to  20  km.  The 
middle  panel  shows,  at  the  entrance  of  the  necking  crust,  part  of  the  energy  is  deflected  from 
its  original  direction  and  no  longer  meets  the  total  reflection  condition.  This  energy  finally 
leaks  into  the  upper  mantle  and  cause  attenuation  of  Lg  energy.  This  can  also  be  seen  from 
the  energy  attenuation  curve  in  the  bottom  panel.  The  model  in  Figure  16  has  a  crustal 
thickening.  Again,  at  the  place  where  the  crust  thickness  decreases  (as  the  wave  leaves  the 
exit  of  the  belly),  part  of  the  energy  is  deflected  from  its  original  direction  and  leaked  into 
the  upper  mantle. 

5.2  Waveguide  models  with  random  structures 

In  addition  to  the  large  scale,  deterministic  structures  discussed  in  the  above  examples, 
small  scale  heterogeneities  also  widely  exist  in  the  crust  and  upper  mantle.  These  small 
scale  structures  can  be  simulated  by  random  velocity  fluctuations  on  a  background  velocity 
and  by  random  interface  undulations.  The  following  models  will  be  used  to  test  how  these 
small  scale  structures  affect  Lg  energy  transport.  Before  applying  the  screen  propagator 
method  to  a  random  velocity  model,  we  first  check  the  validity  of  the  screen  method  by 
comparing  its  result  with  that  from  a  finite-difference  method.  The  upper  panel  of  Figure  17 
gives  a  crust  model  with  5%  RMS  velocity  perturbation.  To  shorten  the  computation  time, 
here  we  use  a  16  km  thin  crust  model.  The  lower  panel  shows  the  comparison  of  relative 
attenuation  curves  between  the  two  methods.  The  solid  line  is  from  a  fourth-order  finite- 
difference  method  and  the  dotted  line  is  from  screen  propagator  method.  The  two  results 
give  reasonably  consistent  results.  This  proves  the  validity  of  the  screen  propagator  applied 
to  the  energy  transfer  and  partition  in  random  crustal  waveguides.  For  this  test  model,  the 
FD  calculation  has  Aa;  =  A z  =  0.125A;ra,  At  =  0.015  sec.  resulting  in  a  CPU  time  of  58 
hours  on  a  SUN  SPARC-4  work  station,  while  the  screen  method  has  Ax  =  Az  =  0.25 km, 
At  =  0.1  sec.,  and  a  CPU  time  0.5  hour  on  the  same  machine.  Both  calculations  have  /0  =  1 
Hz.  For  the  screen  method,  the  cutoff  frequency  fmax  =  2  Hz. 

Shown  in  Figure  18  is  a  normal  crustal  waveguide  similar  to  the  reference  model  except 
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Crustal  model  with  exponential  random  media 
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Comparison  of  energy  attenuation  between  FD(solid  line)  and  GSP(dotted  line) 


Figure  17:  Comparison  of  energy  attenuation  curves  calculated  by  screen  propagator  method 
and  finite-difference  method.  Shown  in  the  top  panel  is  the  velocity  structure  including 
random  heterogeneities  and  in  the  lower  panel  are  the  energy  attenuation  curves. 
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there  is  a  5%  RMS  random  velocity  perturbation  in  the  crust.  The  fluctuation  has  an 
exponential  power  spectrum  with  horizontal  and  vertical  characteristic  scales  of  5.0  and  3.0 
km,  respectively.  Compared  with  Figure  13,  energy  continuously  deflects  from  lower  vertical 
slowness  to  higher  vertical  slowness  due  to  the  scattering  by  small  scale  heterogeneities. 
This  can  be  seen  by  comparing  the  middle  panel  of  Figure  18  with  that  of  Figure  13.  In 
Figure  18  there  is  more  energy  seen  above  the  critical  wavenumber.  The  energy  with  larger 
vertical  slowness  tends  to  leak  into  the  upper  mantle  and  cause  additional  Lg-wave  energy 
attenuation  as  can  be  seen  in  the  lower  panel,  where  the  dotted  line  is  for  the  reference 
model  and  solid  line  is  for  the  random  model.  Figure  19  is  the  energy  distribution  versus 
slowness  and  distance  for  individual  frequency  components.  Comparing  with  those  in  Figure 
14  for  the  reference  model,  we  can  see  that  many  interference  patterns  for  guided  modes  are 
destroyed.  This  is  especially  clear  for  higher  frequency  components.  This  phenomenon  is 
also  related  to  the  characteristic  scale  (or  equivalently,  the  spatial  power  spectrum)  of  the 
random  heterogeneities.  Figure  20  gives  the  attenuation  curves  for  different  characteristic 
scales.  The  upper  panel  is  the  relative  total  energy,  which  is  the  energy  contained  in  the 
whole  seismogram  recorded  on  the  surface  versus  distance.  The  solid  line  is  for  ka  =  1,  the 
dotted  line  is  for  ka  =  10,  and  the  dash  line  is  for  the  reference  model.  We  see  that  for 
the  reference  model,  the  total  energy  remains  basically  constant  beyond  critical  distance, 
which  serves  as  a  checking  point  for  the  numerical  simulations.  The  lower  panel  gives  the 
logarithmic  relative  RMS  Lg  wave  amplitudes,  which  are  calculated  within  the  Lg  window, 
versus  distance.  Again,  the  solid,  dash  and  dotted  lines  are  for  ka  =  1,  ka  =  10  and  the 
reference  model.  In  both  measurements  ka  =  1  lines  give  stronger  attenuation  than  ka  —  10 
lines.  We  see  also  that  the  Lg  amplitudes  attenuate  more  rapidly  than  the  total  energy.  This 
is  due  to  the  effect  of  scattering,  which  diffuses  the  waves  out  of  the  Lg  window. 

The  existence  of  rough  surface  and  interfaces  has  effects  on  both  excitation  and  propa¬ 
gation  of  Lg-waves.  First,  the  rough  surface  near  the  source  region  may  affect  the  source 
energy  partitioning.  Second,  along  the  entire  propagation  path,  the  trapped  energy  in  the 
waveguide  may  be  redirected  to  steeper  angles  due  to  scattering,  causing  it  to  leak  into  the 
upper  mantle.  This  can  cause  additional  Lg-wave  attenuation.  Similar  to  the  example  shown 
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Angular  Spectrum  with  random  media  in  crustal  waveguides  at  f=[0.6,1.9] 
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Figure  18:  Energy  distribution  for  a  crust  model  with  5%  RMS  velocity  perturbations.  From 
top  to  bottom,  crustal  structure;  energy  distribution  versus  vertical  slowness  and  distance; 
energy  attenuation  versus  distance. 
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Figure  20:  Attenuation  curves  for  different  characteristic  scales.  The  upper  panel  is  the 
relative  energy  attenuation,  and  the  lower  panel  is  the  logarithmic  relative  RMS  Lg  wave 
amplitude  attenuation.  The  solid  line  is  for  ka  =  1,  the  dotted  line  is  for  ka  =  10,  and  the 

90 

dash  line  is  for  the  reference  model. 


in  the  last  section,  we  use  a  low  velocity  layer  with  a  rough  lower  interface  to  simulate  the 
effect  of  a  rough  surface.  In  Figure  21,  the  upper  panel  is  the  velocity  structure  including 
a  surface  sedimentary  layer  with  a  random  interface.  The  sedimentarylayer  has  an  average 
thickness  of  1  km  and  a  RMS  depth  perturbation  of  0.2  km  with  correlation  length  of  1  km. 
The  middle  panel  is  the  energy  distribution  versus  distance  and  vertical  slowness.  Compared 
with  Figure  13  for  the  reference  model,  it  clearly  shows,  with  the  existence  of  a  rough  inter¬ 
face,  considerable  energy  moves  from  lower  vertical  slowness  to  higher  vertical  slowness.  In 
other  words,  the  energy  propagation  directions  are  deflected  from  near  horizontal  to  steeper 
directions,  which  makes  more  energy  leak  into  the  upper  mantle  and  causes  extra  Lg  wave 
attenuation.  The  lower  panel  gives  the  energy  attenuation  curves  versus  distance,  in  which 
the  solid  line  is  for  the  model  with  rough  interface  and  the  dotted  line  is  for  the  reference 
model.  It  shows  clearly  the  effect  of  rough  interface  on  Lg  attenuation. 

6  Conclusion 

The  advantages  of  the  half-space  GSP  method  can  be  summarized  as  follows. 

•  Fast  speed:  For  medium  size  2D  Lg  problems,  it  is  2-3  orders  of  magnitude  faster 
than  the  FD  methods.  For  large  distance,  high  frequency  3D  problems,  the  time  saving 
factor  could  be  much  greater. 

•  Memory  saving:  The  GSP  needs  only  to  store  2D  data  arrays  for  each  step  instead 
of  3D  volume  data,  leading  to  huge  memory  savings. 

•  Stability:  The  transversal  Laplacian  is  calculated  by  the  Fourier  method,  therefore 
no  numerical  dispersion  occurs  for  high  frequency  waves. 

•  Intrinsic  attenuation:  Being  a  frequency  domain  method,  it  is  easy  to  incorporate 
various  Q  models  into  the  simulation.  Therefore,  the  method  can  study  the  effects  of 
scattering  and  anelasticity  on  Lg  wave  blockage  and  attenuation. 
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Energy  attenuation  versus  distance  (line-rough  interface) 


Figure  21:  Energy  distribution  for  a  crustal  waveguide  model  with  rough  interface.  From 
top  to  bottom,  crustal  structure;  energy  distribution  versus  vertical  slowness  and  distance; 
and  energy  attenuation  versus  distance. 
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•  Random  heterogeneities  and  rough  interfaces:  Random  heterogeneities  and 
rough  interfaces  can  also  be  easily  incorporated  into  the  model.  The  effects  of  small 
scale  heterogeneities  and  rough  interfaces,  and  their  statistical  characteristics  can  be 
studied  by  numerical  simulations  using  this  method. 

The  use  of  half-space  phase-screen  propagator  has  made  seismogram  synthesis  at  regional 
distances  very  efficient.  Rapid  changes  of  crustal  structure  and  small  scale  heterogeneities 
and  rough  interfaces  in  the  crust  can  be  easily  handled.  The  results  in  this  paper  demon¬ 
strates  the  advantages  and  feasibility  of  the  approach.  P-SV  and  full  3D  problems  will  be 
treated  in  future  work. 
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